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Abstract 

The development of sputtering facilities, at the moment, is mainly pursued through exper- 
imental tests, or simply by expertise in the field, and relies much less on numerical simulation 
of the process environment. This leads to great efforts and empirically roughly optimized so- 
lutions: in fact, the simulation of these devices, at the state of art, is quite good in predicting 
the behavior of single steps of the overall deposition process, but it seems still ahead a full inte- 
gration among the various tools already available for the phenomena involved in a sputter. We 
summarize here the various techniques and codes already available for problems of interest in 
sputtering facilities, and we try to outline the possible features of a comprehensive simulation 
framework, able to integrate the single paradigms in a full simulation, dealing with aspects 
going from the plasma environment up to the distribution and properties of the deposited film, 
not only on the surface of the substrate, but also on the walls of the process chamber. 

Introduction 

The simulation we are going to outline should deal with the specific case of Thin Film Sputter- 
ing processes, within the general framework of Physical Vapour Depositions (PVD). Basically, this 
kind of processes means operating in a ultra High- Vacuum (UHV) environmenio, with tempera- 
tures ranging from RT up to ab. 200°C, occasionally reached in certain parts of the system. A 
preferential direction of investigation will be devoted to Magnetron Sputterers (a strong magnetic 
field - ab. 200 Gauss - is superposed to the static/dynamic electric fields governing the behaviour 
of the plasma) , since the magnetically aided ignition and maintenance of the plasma state has a 
significative improvement compared to non-magnetic technologies, and therefore is currently the 
standard in both research and industrial facilities pi. Main advantages are in fact: 



a few mTorrs of pressure, filled up with neutral species like Ar 
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• reduction of the target-sample distance, and consequently higher efficiency; 

• same ionization efficiencies at lower process gas pressure, which aids both the quality of the 
process and the deposition yield of sputtered atoms; 

• no need to increase the voltage drop, therefore avoiding an increase in the electron mean 
free-path. 

The long-term objective would be the full simulation of a sputter deposition procesfH, which 
means the inclusion of the interdependence, due to their strong interactions, among the different 
phenomena concurring to the global process. In other words, a fully self-consistent, multi-physics 
simulation, where the calculated dynamics of each aspect is integrated with and directly affects all 
the others included in the simulation. A less ambitious approach would, in the first stage, deal with 
the various phenomena separately, neglecting or simplifying their mutual dependence: this would 
look much more like a collection and integration of state-of-art available simulation tools in the 
field (which would become the modules of the global approach), with the final scope of merging 
them into a framework, able to efficiently provide parameters, derived from one simulation module, 
as inputs to the modules directly related to it. 

A straightforward novel application would be studying and predicting the interference of side- 
apparatus with the core sputtering facility. Specifically, this general approach could be applied to 
shields, simple or engineered ones, invasive sensors, non-standard materials employed inside a sput- 
tering facility. These insertions rely currently more on experimental trials, rather than numerical 
simulation, which leads to empirically, roughly optimized solutions (given the extreme heterogeneity 
of the facilities available in this sector) and the need for new experimental tests - whenever changes 
to the facilities are made. 

Expectations from the results of the simulation are more a hint about how the facility could react 
upon the insertions of non-standard features, and certainly not an ab-initio, stand-alone calculation 
of optimizing parameters, given the assumption that 

treatment of a real plasma discharge in a reactor (is) not easily amenable to analytical 
or even numerical approaches, due to the variety of complex interactive phenomena f4l 

Therefore, an interesting approach would be the progressive experimental validation of results from 
the simulation, and re-tuning of the parameters involved, until obtaining a good agreement with 
well-known, or easy-to-check, experimental data. Once this first result is established, one could 
move further to use simulation data as a realistic base for the optimization of currently available - 
as well as upcoming - technologies. 



1 The Plasma source 

As stated before, this is probably the most difficult part to simulatclfl. In fact it involves the dy- 
namic simulation of motion and collisions (with appropriate cross-sections) between plasma particles 
(heavy ions and electrons) and background gas neutral atoms, with a particular focus on elastic and 

first attempt in this direction, with a success which is up to now very Umited, has already been done with the 
development of the NEPTUNE Sputter code [T] [23] 

•^for a more general review of the approaches to the problem, which we synthesize here, we suggest I16| 



2 



ionization collisions in the plasma environment. Parameters describing the plasmas can be highly 
variable, but gas pressures around 1 Pa, with electron densities in the order of 10^" cm~^, should 
be considered plausible. This requires an algorithm including a self-consistent simulation of the 
plasma charged particles motion, due to: 

• electric fields (mainly dc/rf, depending on the facility, with potential of 200 — 400 V between 
cathode and anode; plus turbulent electric fields occurring in the electron trap region [3], 
whose strength is a few V/cm) and 

• magnetic fields (in the case of Magnetrons: these can always be considered static, with an 
intensity of about = 200 Gauss), 

including the mutual interactions between the sources of these two. Given the physics introduced 
above, it is a very used paradigm - in the general description of the plasma [8] - the introduction 
of continuity and Poisson equations: 



dt 



V ■ Jj(e) — Ri{e) (1) 



V2(/) = -^(n, -ne) (2) 
Co 

to match the behaviour of the ions (electrons) - in eqU] - each with its own density njfg-) , with 
the determination of electrical potential - in eqj^l Notice that, given the ionization processes 
occurring in the plasmsl^, the continuity eqs. must embed a source term R. 

The plasma environment has normally temperatures as low as possible, since hotter plasmas are 
more difficult to control, however, peak temperatures of the residual gas can reach also 900 K or 
more [261. 



1.1 The external fields 

A first step in the modeling of the overall 'plasma ignition' process relies upon an accurate cal- 
culation of the magnetic field generated by the permanent magnets inside the chamber and, in 
particular, in the proximity of the target's surface, exposed to the action of the plasma (the so 
called ^sheath' region). The flux density of the magnets is a known parameter, as well as their 
position, normally a few cm behind the cathode-target. 

It has been found in [9] that the particles' trajectories mostly affected by the action of the magnetic 
field are indeed the electrons' trajectories E|. This is why the accuracy of the field calculation must 
be highei0 where the density of the plasma electrons is maximal. The solution of this magnetostatic 
problem should pose no particular difficulty, and can be afforded i.e. with FEMLAB modeling, as 
it was done in [25] . 

An alternative approach is to introduce an effective potential (^) picture: 

^^ ^^%-^r^v # (3) 



^Which will be described in detail in the next paragraph 
^since ions in the plasma are not significantly magnetized 

^this would mean a finer grid in the zone up to ab. 4 cm away from the cathode, or even the calculation to be 
performed in this zone only 
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where q and m are the electrons' charge and mass, Pe its canonical momentum, whose position is 
described in a cylindrical systeirQ of coordinates (r, z,9), under the action of an electric potential 
(j) and a magnetic vector potential Aq. In particular, the z direction is the one along the target- 
substrate distance. In this way, one can make a rough description of the zones where most electrons 
are trapped, rather than, calculating the position and extension of the sheath and pre-sheatJ^ 
regions, using approximations: 

• for Ag, given basic characteristics of the facility simulated, like the position of the target and 
the applied external magnetic field; 

• for (f) in the sheath, instead, referring to theoretical models like Child's law [TT| : 

• and finally (j) in the pre-sheath is normally considered uniform [9 , or the smalU potential drop 
in this zone can be accurately measured [12j and is in the order of a few V. 

Once ^ is known, the next step is the application of a Hamiltonian approach, like the one suggested 
in pi51, to derive the dynamic properties of the particles to be followed. 

As a first approach to validate the simulating code, the surfaces of the chamber can be considered 
as perfectly conductive, grounded sheets, so to deal with an easier situation. A proposal to draw a 
more realistic picture in case of non-conductive insertions will be explained in the following. 
In conclusion, the intensities of these external fields should be treated as known conditions (i.e. input 
parameters) in the simulation, since the sources are known and controlled within a high degree of 
approximation. What is left to calculate (or directly measure) is the intensity of these fields within 
the real geometry of the equipment: the chamber where the process occurs, plus eventual insertions. 



1.2 Particles' Trajectories 

As emphasized in [7 , the general problem of particles' trajectories in a sputtering simulation can 
be split in at least two different frameworks, which will be treated in an almost independent way: 
electrons and ions. Neutrals will be included as static targets of the collisions occurring in the 
plasma, rather than as dynamic particles themselves. 

Electrons of interest in the plasma dynamics are the so called fast electron^^ These are mostly 
generated in the sheath region, or in the proximity of the cathode [3], by the combined action of 
the external fields and of the collisions with neutral particles. Some other are also the secondary 
electrons (also known as bulk electrons or 'discharge electrons') generated as a side effect of the 
impact upon the target of the heavy positive ions (or to a minor extent, within the interaction with 
neutrals). It is possible to describe the phenomenon by an empirical yield 7, so that the secondary 
electrons current density je at the target (position on the z-axis): 

Je(0) = -7J.(0) (4) 

is related to the ions' current density ji. 

The motion of the electrons can be thus described starting from the integration of the equations of 
motion: 

i^ — (E + vxB] (5) 
m \ / 

^particularly indicated for the case of a planar Magnetron 

*the zone immediately further away from the cathode, where the tail of most energetic ions is generated 
^If compared to the one occurring in the sheath region 
^"i.e. those sufficiently energetic to ionize neutral atoms of the gas 
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Once an efficient integrator for cq.(l5]) has been provided, one is able to track with high accuracy 
the location (i.e. orbits), velocity and kinetic energy of the electrons. In this way, one can improve 
the qualitative picture, given by the potential 'I', analyzing as the first this collisionless motion of 
the single electrons in the plasma!"]. 

Naturally, in order to derive realistic results, it is necessary to include the collisions among the 
electrons and the other species inside the sputtering facility, since the internal pressure is not so 
low to neglect them as in thermal evaporators. This will be the topic of the next subsection. 

In a sputter, ions are originated by some of the collisions occurring in the plasma^, which 
means that these sites can be accurately predicted !7 , and lay in the electron trap region. Once the 
ionization has occurred, the calculation of ions' trajectories follows the same guidelines introduced 
for the electrons' case. 

This means that the ions are subjected mainly to external fields (as already outlined, in particular 
to the dc/rf field, since the magnetic field has a reduced influence on ions, and could also be 
neglected). Apart from this major contribution in the r and z directions, experimental data |12| . 
have emphasized a noticeable amount of random energy in the direction. For this additional 
feature, normally turbulent electric fieldJ^. or further collisions occurring in the plasma, are held 
responsible. 

Finally, notice that, within the drift-diffusion framework leading to eqUJ one can also describe 
globally the motion of both positive/negative charged particles in terms of a flux J [5], which is, 
including the action of the magnetic fleld for the specific case of a DC Magnetron: 



Ja = {-ha yj^// + fl^E^ + X h) 



(6) 



where h is the unit vector in B direction, and for each charged particle (a := {j,e} ), we have 
defined its average rate of collisions with neutrals (n) , function of the density of the residual gas 
n, and therefore its effective mass m* , the components of the mobility: 



Ma 



Ma 



B 



Mi 



M 



// 



(7) 



and of the diffusion coefficients: 
£)// ^ ^-sTg . 



d" 



(8) 



for a certain cyclotron frequency Wa = laB/nia. 

^'^The calculation, for the case of a planar Magnetron, will indeed show that some of these orbits are confined near 
the cathode by the shape of while others are unconfined and the corresponding electrons are able to 'escape' in 
direction of the anode, before releasing most of their energy in collisions. This collisionless description in Plasma 
Physics is a well established formalism also known as 'Vlasov-Maxwell equations' 

^■^For further details see next subsection 

^^even if they are about 4 times smaller than the dc field already in the prc-shcath region 
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1.3 Collision phenomena 



A variety of different kinds of collisions can occur in the plasma environment. A rough classification 
is given in the following. 

1. Electrons VS Neutrals. They can be of three main typologies: ionization, elastic scattering, 
excitation. The firsts are probably the most fundamental, since they generate and maintain 
the plasma itself: every time one occurs, the corresponding electron loses an energy equal 
to the ionization potential of the neutral specie, plus the kinetic energy of the (secondary) 
electron emitted. Elastic scattering produces a reduced energy loss of: 

= sin^ia 2) 9) 

Finally, there is the contribution from the excitation collisions. These last ones have been 
shown [13] to be ab. 1% of the total collisions, so are often neglected, along with Penning or 
double ionization phenomena. 

It must be remembered that each kind of collision has its specific scattering cross section 
da/dTl, which varies with the energy of the incident electron. Detailed tabulated values of 
all three cross-sections are available in the literature [5 and can be used in the simulation; 
interestingly, for highly energetic impactj^. these cross sections are peaked at small angles 
whatever their type is. This means, most of the collisions occur without a substantial change 
in the canonical momentum of the electrorF^. 

2. Electrons VS Charged Particles (or: Coulomb collisions). Considering a typical Magnetron 
discharge, the frequency for these scattering events can be evaluated about five orders of 
magnitude less than the previous case [7]. This is why most of the models for the sputtering 
behaviour neglect this kind of collisions. 

3. Ions VS Neutrals. This kind of collisions alter significantly the trajectory of the involved ions, 
and therefore are important to understand the 'landing sites' of the ions on the target, and 
consequently the scattering process. In fact, these collisions (both elastic and charge exchange 
ones) are held co-responsible for the azimuthal randomness introduced before. To understand 
the role of these collisions, it has been estimated that an average ion created in the pre-sheath 
region of a Magnetron sputter has about 35% probability to collide elastically with a neutral 
particle before reaching the target [9]. 

4. Ions VS Ions. These phenomena, which in principle occur with a significant probability, can 
be neglected for energetic reasons. In fact, in this process the energy transfer rate is much 
slower than the ion loss ratcl^. 

For this kind of simulations, Monte Carlo (MC) approaches are considered the standard in state-of- 
art literature, with a considerable variety of particular cases. Generally speaking, MC codes used 
in this case must be able to follow a large number of particles at a timc£l, keeping track of their 

^''Which means, with kinetic energy X > 60 eV, which often occurs in a typical Magnetron 

^^And this can be interpreted as a long permanence time for the electrons, compared to the average time between 
collisions, before they are scattered away from the trap 

^®As it can be verified from typical parameters of sputtering facilities, see i.e. [4] 
remember the densities are of about 10^" particles/cm'^ within a cylindrical zone, whose typical dimensions are 
5 cm (for both radius and height) 
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cinematic parameters and providing an accurate mechanism to simulate the recurrence of colhsions, 
each with the adequate cross section. UsuaUy, initial conditions for the particles should be unbiased 
in the code. 

The first and rougher approach is to make use of the approximations cited above, plus an intuitive 
reduction of the number of particles involved in the simulation. This strategy relies on the con- 
sideration that, since charge/mass ratio remains the same, essential physics can be captured with 
a much smaller number of particles than those in a real plasma |16| . therefore, some simulations 
have been run with just 10'^ charged particles. 

The usual way to introduce collisions in this numerical scheme is to divide the motion of the par- 
ticles in time steps, and for each generate random numbers to be compared with the probability of 
a collision to occur, referring to tabulated values for the total cross sections available, for example, 
in [17] . One afterwards discriminates between the various kinds of collisions by comparison among 
the different cross sections related to each case, again by generation of random numbers to decide 
whether the collision was i.e. elastic or ionizing. After that, the energy of the particles involved 
is varied in accordance to the rules stated above (and in [7]), and the motion after the scattering 
is altered in correspondence. Approaches of this kind has led to a noticeable variety of slightly 
different models, which go under the name of "Monte Carlo collision" (MCC) algorithms 
A somehow more intelligent use of computer resources is based upon the so called " Particle- in-cell" 
(PIC) modeling. This type of alternative, competing algorithms deals with the problem in the 6D 
parameter space of the particles, subdivided by a grid in a large set of cells, which are 'efficiently 
sampled' by the particles themselves [B]. This means that we are making use of two main ideas: 

• if we ideally divide the phase space in cells, most of them are unoccupied most of the time, 
so that it is useless to include them in the calculation, complicating the resolution even in a 
Vlasov-Maxwell approach: therefore, each step of the simulation brings along a reduction in 
the number of cells which will be effectively treated in the algorithm at the next step; 

• to perform a full simulation of the interactions among particles, it is not needed to calculate 
them all directly; rather than, it is possible to follow an averaged approach in 3 calculation 
steps for each time step: i) once the coordinates of the particles are known, they are used to 
calculate currents and charges deriving from the particles themselves; ii) these equations are 
integrated on the grid, so to obtain values for the electromagnetic fields in each cell, which 
iii) will be used to extract the forces used in the next time step of the equations of motion 
integration. 

These two ideas, though quite intuitive, have proven highly effective in reducing the computational 
time required for a full simulatiorF^. It is useful to remind that PIC algorithms are not necessarily 
collisionless, as they can introduce collision phenomena by the 'finite-size particles' scheme [TB] . 
A further improvement in the PIC approach is the so called 'implicit' version. In fact, within 
the traditional^ PIC algorithm described above, field equations need only the sources from the 
previous time cycle - and the equations of motions need only the fields from the previous time cycle 
- a simple procedure, whose limit is the severe stability conditions to be matched (see i.e. [18|). 
The power of the implicit approach, instead, relies in the implicit coupling, introduced between the 

'^^In example, it can be calculated that replacing the full naive 3D approach for 10* particles with a 64x64x64 cells 
simulation was able to reduce the computation time on a 10 Tflops cluster of about 6 orders of magnitude 
^^Which in turn can be named the 'explicit' version 
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discretization of the equations of motion for the particlej^ 

ri+l 

D 



— V + 
P P ^ 



(e;'+''(x^+i/2) + X b;j(x^+i/2)) (10) 



nip 

and the discretization of the Maxweh equations for the fields: 

VxE"+« = -l?— 
c At 

c At c 

V . E"+^ = -47rp"+^ 

V • B"+^ = (11) 

Notice that the intermediate values for whatever variable, say Q, are computed as: 
Qn+9 = (1 _ 0)Qn 6ig"+\ whcrc 6 € [1/2, 1], since for 9 < 1/2 the algorithm is known to be 
unstable [19j . The main advantage deriving from these coupled equationto is in the reduction to 
only a set of coupled fluid moment and field equations to be solved. Within a single time-step 
there is no more need of iterations to compute all of the parameters in the model. For a thorough 
description of the approximations introduced, to provide accurate solutions for the coupling, we 
refer to P ] and [TO ] . 

We emphasize further that a successful implementation of these PIC schemes already exists and is 
provided by the CELESTE3D© codcF^. of which several releases are available. 

An interesting consideration to understand, for all the aspects of the plasma source described 
above, is: whether conditions are not too extreme to be treated within a fluid model, is there a 
chance to adopt multi-physics software^? And by comparison of the results obtained by the soft- 
ware, with those from state-of-art MC techniques, does it exist a valid range of conditions, where 
such softwares can provide useful results with a considerable reduction of computation time and/or 
resources required? 

The answer, up to the studies already available in the literature, can be considered carefully posi- 
tive. For sure, in fact, a successful attempt in modeling (by a hybrid approach partially based on 
COMSOL modules) several parameters of a plasma discharge in a Magnetron device has been made 
in [26], obtaining a good match with available experimental results. The assumptions made in the 
cited model are essentially: 

• the sources (ionization sites) in the plasma and the sampled particles' trajectories are calcu- 
lated by a MC approach like the one suggested in [S], under the same assumptions; 

• gas heating effect^^ are taken into account in the simulation, and lead to rarefaction of the 
gas itself, which is modeled as ideal; 



^"whcre p is the particle index, 9 is the discretized time-step, n is the iteration step 

^^in fact, the presence of Ep+® in eg.s IIUI requires the knowledge of the advanced electric field to calculate the 
particles' position, and viceversa, since the charges/densities rely on particles' sources in llll 
http: / / code.google.com /p/ celeste / 

Such as i.e. the COMSOL® ' Plasma module' (hereafter the ©will be omitted for brevity), eventually combined 
with other modules 

^*due to the energy lost by charged particles collisions with neutral particles of the residual gas 
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• data from the two points above are used as input parameters for the fluid-Poisson model, 
implemented in a COMSOL Multiphysics environment, solving the continuity and Poisson 
equations; 

• the full 3D problem has been reduced, given symmetry propertiej^. to a 2D mesh perpen- 
dicular to the cathode surface, and properly discretized. 

Additionally, the boundary conditions set for the equations ([1]) and ^ are: 

(for the electrons) := (at cathode and walls) 

(for the ions) rii := (at walls, or defined by) eq. (|1]) (at cathode) 

(j) := (at walls, grounded and) — —Vappi (at cathode) (12) 

Notice that actually the only bound self-consistent with the MC calculation is for the cathodic ions' 
density: all the others rely on the assumption of a perfectly conducting, grounded surface of the 
walls. 

The problem is afforded by the use of stationary 'Direct UMFPACK' solverQ, both linear and 
non-linear, to calculate the magnetic field inside the process chamber (with Magnetostatic modules) 
and the COMSOL Chemical Engineering Module for the solution of eqs. ([1]) and ([2]). The authors 
claim the results were in good agreement with experimentally available data from Langmuir probes, 
concerning the temperature of the residual gas, the electrons' densities and ions' flux, and finally 
the electrostatic potential in the chamber. 

Results from the work above, by the way, should be considered carefully. A first major issue is that, 
in fact, the high ratio of the parallel magnetic field to cross field to electron mobility in a Magnetron 
is very high, and the fluid model involved in the Plasma module code barely handles such extreme 
conditions. Such a module could thus be applied to the case of non-magnetically-aided facilities 
only, whereas for Magnetron-like environments hybrid approaches must be invoked. Another point 
to emphasize is the possible difficulty in dealing with RF sputtering (which up to now has been 
recalled quite a few) , especially for the case of recent RF sputtering facilities employing Inductively 
Coupled Plasmas (ICP), because of a series of technical advantage^. 

ICP sputtering is simpler to simulate from some points of view: i) at pressures much lower than 
1 mTorr, transport of ions can be considered collisionless, instead of diffusive; ii) because of the 
low plasma potential, ion bombardment of wall surfaces is reduced, and therefore the possibility of 
co-sputtering and charging effectf[^ is less significative. In any case, even if these aspects require 
less accuracy in the simulation, it turns out that the behaviour of the plasma itself can be tricky to 
simulate. In fact, in literature are reported cases of failure in the predictions about basic parameters 
of the plasma (like the electron density), based upon simulations run with the 'Plasma' and the 
'RF' modules from COMSOL. A specific example is given by the low pressure regime of a jet ICP 
reactor, where important discrepancies were found between experimental results and the numerical 

^^recalling the azimuthal non uniformity introduced in the previous paragraph, this should be considered a good 
approximation more than an exact hypothesis 

http: / / www. cise. ufi. edu / research /sparse/ umfpack / 

^'^I.e. they avoid contact of metaUic electrodes with the plasma and ionization mechanisms on isolating or con- 
ducting walls (the plasma is heated by the RF field induced by an external antenna), which leads to a high degree 
of plasma purity, lower pressures are possible, making a more directional technique available, and finally the absence 
of high voltage sheaths avoids potential co-sputtering from the walls. 

^*In fact, some ICP sputterers use even dielectric shields [4]; a more detailed discussion of these phenomena will 
be presented below 
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simulation obtained for the plasma sources, based upon interaction between the electromagnetic 
fields and the process-gas flow [27] . 



2 The sputtering process 

This process deals with the dynamics of energetic ions, which are accelerated towards the plane 
surface of the target material. The impact with the target crystal generates recoil cascades inside 
the lattice, displacing a certain number of atoms from their equilibrium sites. Those atoms on the 
target surface, which are not confined along the external direction, are hold on their sites from the 
Surface Binding Energy (Ef,), which is usually lower than the Displacement Energy of the bulk 
material [3]. Whenever these surface atoms acquire, from the cascades, an energy higher than Ei,, 
they are extracted from the surface itself [2], and this leads to their sputtering away from the target. 
Because of the physics of the process, a thin laye^ of the target material is enough to perform 
an accurate simulation, reducing the time spent calculating cascades which will not contribute to 
sputtering. The most important value in describing the result of the sputtering process is the 
sputtering yield Y, simply defined as: 

sputtered atoms 

Y = ■■ 13 

incident ions 

where clearly, the higher the yield, the most effective is the ionic 'bombardment' for the deposition 
purposes. The yield itself can be obtained experimentally, or derived from semi-empirical models, 
like the Sigmund, Bohdanski, Yamamura or Wilhelm formulas Basically, there are two different 
approaches in solving the problem of simulating the sputtering yield and the energetic/angular 
distribution of the sputtered atoms. 

A first choice is relying on experimental data, which are accurate and available for a consider- 
able amount of different materials, in order to opportunely tune the parameters of semi- analytical 
models. In particular, the probability for a scattered atom to leave the surface of the target has a 
distribution dJ/dE with the shape [H) : 

dE - "^WTW ^''^ 

rather than, as an input to the 'transport module' (see the next paragraph), it is possible to express 
the starting energy Eq of each simulated atom as [12] : 

f E, 

= TJ2 (15) 

where Ei,om is the energy of the incident iorFl. 

while is a random number indicating the probability for an atom with energy Eq to leave the 
target. A similar approach is used for the angular distributions, which are treated as respecting the 



^^Usually within the range 30-100 nm, from heavier to hghter ions, since these last ones tend to penetrate the 
target in deep |15l 

^''which can be derived from further experimental data or as an input of the 'plasma source' environment 
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Knudsen cosine law {dJ/d9 — Dcos9), for the starting zenith angle 9o [35], and full randomness, 
for the azimuthal starting angle: 



where and are other random numbers. The starting trajectories sampled out in this way can 
be considered a rough approximation of the true input for the transport module, leading in any 
case to good predictions 22 . 

A more refined version of this approach involves instead the full treatment of the collisional 
dynamics in the solid, in order to retrieve the parameters above j2l], [S^. The standard is here 
again the adoption of MC algorithms, and in particular the TRIM (TRansport of fons in Matter) 
module of the SRIM packagj^. The core of the program is a detailed tracking of the incident 
ions and recoil processes, until their energy falls below a threshold indicating they are not anymore 
candidates for sputtering contributions ,15i . Given that no assumptions are made about the lattice 
structure of the target, the code is reliable also in the case of amorphous material^. Statc-of- 
art implementation of the code makes use of the simpler but robust 'hard sphere' approximations 
feallSp for the treatment of low-energy recoil events; more accurate theoretical tools are indeed 
available (like the 'Universal potential' of eqUTJ the 'Kr-C or the recent 'ZBL electronic stopping' 
approaches [H]), and the last one in particular is used - in the TRIM package itself - for high energy 
events. 

Basically, the TRIM code requires four input values for the simulation to be performecj^: 

• the lattice displacement energy and binding energy; 

• the surface binding energy; 

• the sputtering yield for normal incidence. 

Moreover, the simulation allows to set a certain number of impacting ions, with a specific energetic 
and angular distribution: the effect of the angular incidence on the corresponding normal sputtering 
yield will be considered. Given these data, TRIM is able to calculate as output the sputtering yield 
dependency on ions' energy, the number of atoms displaced but not able to reach the surface of the 
target, and the angular distribution of the sputtered atoms [TS]. However, it must be remembered 
that, if energetic plots have since a long time proved to provide good predictions - once accurately 
calibrated [21] - the angular plots (and therefore the direction of emitted atoms) seem to suffer from 
some intrinsic limit of the code [24]. A possible hint, about the reason for this discrepancy with 
experimental data, could be the roughening of the surfacj^. which both increases the yield - and 
this can be taken in account by minor adjustments - but significantly alters the angular distribution 
too, in a much less predictable way. 

An alternative to the TRIM code is the implementation of dynamic MC modules (i.e. the so 
called SASA-sp: the sputtering version of the dynamic-SASAMAL code), which are used in the 

^'^further details, papers and releases available under http://www.srim.org/ 

rather than, polymeric materials of interest to us should not pose issues about the structure of the algorithm 
itself 

Actually the code is provided with a library including over 28,000 different cases. However, since the package 
was calibrated for a specific case |24) . significant discrepancies with experimental results could be obtained in case 
no further calibration is performed 

^"^Which is ignored from the TRIM code, and increases as the sputtering process continues. 



6o = arcsin^e 
4>o = 27r^0 



(16) 
(17) 
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'NEPTUNE-Sputter' simulation environment [23]. In this case, the incident ion trajectory is seen 
as a set of straight-hne segments of a length, which depends on the local atom density of the target 
material. This density is dynamically varied as the process goes on. At every vertex, a collision 
takes place with a certain specie of atom, according to the relative concentrations and scattering 
cross sections in the zone where the collision occurs. The electronic energy loss is subtracted from 
the kinetic energy of the incident ion, and a new motion/collision step in the algorithm is performed 
[23j . However, even if this code should have the enormous advantage of being part of an almost full 
simulation of sputtering facilities, it must be reported that it has been rarely cited - and therefore 
tested - up to now, in available literatur^fl; among the few examples, see [T] (chap. 3). 

Beside the traditional interest in the sputtering simulation from the target, an under- investigated 
problem is the role of the ions' impact against the walls or other side-apparatus (sometimes made 
of non-standard materials) inside the deposition chamber, which i.e. could lead to co-sputtering 
phenomena. In the specific case of amorphous dielectric structures used for shielding 2 , charging 
effects of the shields could also be possible. If so, it would be very interesting to understand the 
contribution from this charging to the global electromagnetic environment where the process occurs: 
this interdependence could make multi-physics approaches particularly adapt for these studies El. 

3 The transport process 

The process of transport of the ejected particles, through the low pressure process gas, makes use of 
much of the formalism introduced for the plasma environment. An additional simplification allows 
to ignore the interaction with the applied electromagnetic fields: in very good approximatiorF^ the 
sputtered atoms are all neutral [3j. 

Therefore, this step of the simulation essentially deals with the interaction among sputtered atoms, 
plasma particles and the background gas atoms. This scattering in gas phase can be modeled 
through different approaches, which can be roughly listed as follows, in order of increasing refine- 
ment: 

1. hard-sphere approximation, which essentially relies on the simple definition of an average 
'radius of interaction' r^, leading to a step potential: 



this is especially used to model elastic collisions; 
2. Born-Mayer interatomic potentials, which provide a range of interaction among the particles: 



''^This could be due to the absence of an open-source code, which is the case for the TRIM code, or the reduced 
library of cases-of-study and optimizations 

^''Considering in example the COMSOL framework, it is not clear to us whereas any of its modules could be able 
to directly deal with the sputtering process problem. In fact, up to the knowledge of the author, applications of the 
COMSOL package in atomic extraction from surfaces have currently dealt with evaporative sources ([28], or |25) . 
where the COMSOL simulation is indeed limited to heat transfer effects) or CVD sources only (where the nature of 
the process is mainly chemical). 

^'^For a typical clean metal or semiconductor surface, the percentage of charged emitted particles is only Ri 10"'* 
of the total 




if |ri - r2| > r„ 
Vint if |ri - r2| <ro- 



(18) 




(19) 
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a thorough investigation of plausible values of the parameters A, B, has been given in 

HQl; 



3. Lennard-Jones potentials, of typical use in condensed matter physics, with the form: 



V = AXl, 



r ' \r 



(20) 



4. a modified version |21] of the so called Universal interatomic potential, which has the forr 

y = -^>^L4-)'+i'-^)xuir) (21) 

m 



r J \ r 



where the 'universal screening function' xu has been usei 
5. the Firsov interatomic potential, again based on a screen functiorF^ V'('') 



V={^^)i;{r/x) (22) 



Whatever the particular inter-atomic potential chosen for the simulation, the scattering angle f3 
can be classically calculated from the integral definition: 

P = -2vrp /^°° '^''^''^ (23) 

ir„,„ - Vir)/EcM - pVr^ 

where p is the collision parameter. The calculation of the integral ([25)1 reveals some issues which 
can be overcome by the application of proper strategies [21] . 

Once accurate potentials have been introduced to describe the interactions of the scattered atoms 
with the residual gas, and the relative cross sections are known from tabulated or experimental 
data (see par. II. 3|) . the simulation will be able to take into account both the thermalization and 
diffusion processes, along with the direct fiow, involved in the transport of the sputtered species. 
A typical simulation scheme is the following. MC algorithms are used to introduce and calculate the 
effects (as described above) of scattering events in the direct fiow of the emitted particles, at a rate 
dependent on the mean free path of the particles at a certain pressure, as available from experimental 
data [22]. Once a particle has been thermalized by the collisions (e.g. its energy is within a a from 
the thermal energy ksT), one can decide to keep tracking the particle by MC techniques, or rather 
spare computational resources switching to analytic diffusion equations. Interactions with charged 
particles and inelastic scattering are usually neg 

lectetS- 

The implementation of this transport step 
alone, given plausible description of the starting locations for the atoms, has proved able to provide 



^^the van der Waals attractive contribution is introduced to fit the behavior at impact energies lower than 2 eV 
■^^it is semi-empirically defined as 

Xu = 0.1818e-3-2^ + 0.5099e-f'-S''23x _,_ o.2802e-o *''28^ + 0.02817e-0-20i6^ 

where the reduced distance is a; = r(Zj'-2^ + Z^-^^)/0.89ct 

^"the so called Nikulin screen function, for more details see [22) 

'^-^ Because ion and electron concentrations are small in comparison with the atom quantity (the range of ionization 
is 10~^), and average sputter atom energy (about 10 eV) is not enough to excite and ionize background gas atoms [3] 
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results in good agreement with experimental data [H] [52], especially concerning the thickness, 
angular and energetic distribution of the sputtered layer in the zone of the substrate. With the 
adoption of the most refined interactions, moreover, the accuracy was fully satisfying even for the 
case of multi-target deposition^. It is not fully clear, though, if the calculation of thicknesses on 
the chamber walls (e.g. far away from the substrate), explicitly shown in some results j21j . can be 
considered reliable or not. 

An interesting point yet to investigate is whether a multiphysics approaclF^ would be able to embed 
one or several of the approaches above, with results at least comparable to those obtained from MC 
calculations, once the appropriate collision model has been defined in the software, . 



4 The deposition (Film growth) 

The effort in simulating the various steps of a sputtering process has, up to now, essentially ignored 
a properly contextual description of the deposition. In fact, since the transport of the sputtered 
species is essentially a physical process, and most of them are furthermore thermalized already 
before reaching the substrate |21| . no particular care about the surface chemistry is required, and 
the angular distribution is considered enough to understand how uniform the final film will be. 
Another point of view is the thorough analysis of mechanical/electric/thermal properties of the 
deposited film itself |3j, and therefore of the (micro)structures obtained on the substrate, neglecting 
any reference to the specific parameters of the deposition process. 

However, we emphasize that it would be interesting to understand, even at a low-detail level, 
the interaction among the neutral transported species, on one side, and the substrate material 
(semiconductors, metals, plastics), where the accommodation on the final surface occurs. The main 
aim is to have a clue of the mechanical stress of the deposited material, depending on its thickness 
and deposition-rate, upon surfaces of different geometries or chemical-physical properties. An 
immediate, important application would be the optimization of the geometry and/or the material 
used for shielding solutions, in order to reduce the issue of flaking phenomena [5]. For these 
purposes, it is envisaged that a multi-physics approach could be ideal in order to conjugate the 
effects of mechanical stress (deriving from the increasing, eventually non-uniform thickness of the 
deposited film), and the thermal stress (induced by the relatively hot plasma environment - by 
irradiation - and by the kinetic energy of the atoms impacting on the substrate surface). Given 
the nature of the phenomena to include in this last section, it should be considered the aspect 
of the global process less necessary to be integrated with the others. We claim, in fact, that an 
independent module specifically dealing with the film properties, according to data retrieved from 
the other three modules, would already be very useful to cope with this last aspect. 



Conclusions 

In this work it has been given a thorough review of previously made attempts for the simulation 
of single phenomena in the field of sputtering processes. For each case, the methods available 
have been discussed, emphasizing their limits and perspectives. Whenever available, considerations 
about the effectiveness and the agreement with experimental data have been reported. 
In addition, a complete framework embedding the various phenomena involved in a sputtering 

^■^With remarkably good results also for the final stoichiometry of the film |22| 

*^I.e. a 'Particle Tracing' module, example again taken from the COMSOL software framework 
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facility has been given, and the problem for a comprehensive interaction among available procedures 
is posed. Where plausible, adoption of features available in specific multi-physics software has been 
suggested as a possible strategy to simplify the simulation, compared to state-of-art numerical 
approaches. Given the current fragmentation of the approaches dedicated to the various aspects of 
the sputtering, this alternative strategy is considered worth to be investigated further. 
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